function [v s theta]=euler_convert(lat,lon,w); 


%function v=convert(lat,lon,w); 
%
%Converts a 3d rotation vector to a velocity at a point on Earth.
%
%INPUT
%
% lat, lon, of a point on Earth (degrees)
% w,        Euler rotation pole in lat (degrees), lon (degrees), and angular velocity (degrees per million years)
%
%OUTPUT
%
% v,     velocity vector (mm/yr) oriented along North, East, and downward dimensions. 
% s,     speed (mm/yr)
% theta, direction in degrees counterclockwise from North. 
%
%After Cox and Hart, 1986.

if nargin<3,
  help euler_convert; 
  return;
end;

%convert location of moving point to a cartesian unit vector
lambda=lat*pi/180;  %latitude,  in radians
phi=lon*pi/180;   %longitude, in radians
p(1)=cos(lambda)*cos(phi); %x-direction
p(2)=cos(lambda)*sin(phi); %y-direciton
p(3)=sin(lambda);

%radius of small circle
ro=6378e6; %radius in mm 
           %r=sind(distance(lat,lon,w(1),w(2)))*ro;

%convert euler pole to cartesion in mm/year
w(1)=w(1)*pi/180;
w(2)=w(2)*pi/180;
w(3)=w(3)*1e-06*pi/180; 
E=[cos(w(1))*cos(w(2)) cos(w(1))*sin(w(2)) sin(w(1))]*w(3);

%calculate cross product
vc=ro*cross(E,p);

%convert from cartesian rotation vector to a local velocity
T=[-sin(lambda)*cos(phi) -sin(lambda)*sin(phi) cos(lambda);
   -sin(phi)              cos(phi)             0;
   -cos(lambda)*cos(phi) -cos(lambda)*sin(phi) -sin(lambda)];
v=T*vc'; 

%speed and direction
if nargout>1,
  s=sqrt(sum(v.^2));
  theta=atan2(v(2),v(1));
end;
